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Abstract. The profiles of narrow lattice solitons are calculated analytically using 
perturbation analysis. A stability analysis shows that solitons centered at a lattice 
(potential) maximum or saddle point are unstable, as they drift toward the nearest 
lattice minimum. This instability can, however, be so weak that the soliton is 
"mathematically unstable" but "physically stable" . Stability of solitons centered at a 
lattice minimum depends on the dimension of the problem and on the nonlinearity. In 
the subcritical and supercritical cases, the lattice does not affect the stability, leaving 
the solitons stable and unstable, respectively. In contrast, in the critical case (e.g., a 
cubic nonlinearity in two transverse dimensions), the lattice stabilizes the (previously 
unstable) solitons. The stability in this case can be so weak, however, that the soliton 
is "mathematically stable" but "physically unstable" . 
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1. Introduction 

Solitons are localized waves that propagate in nonlinear media where dispersion and/or 
diffraction are present. They appear in various fields of physics such as nonlinear optics, 
Bose-Einstein Condensates (BEC), plasma physics, solid state physics and water waves. 
The dynamics of solitons is modeled by the Nonlinear Schrodinger equation (NLS) in the 
context of nonlinear optics which is also known as the Gross-Pit aevskii (GP) equation 
in the context of BEC. 

In the study of stability of solitons in a homogeneous medium, it is useful to consider 
the d-dimensional focusing NLS 

iA g (z,x) = -VA-IA^A, (1) 

where z is the longitudinal coordinate, x the coordinates in the 

transverse plane, V 2 = d XlXl + • ■ ■ + d XdXd is the Laplacian operator and the nonlinearity 
is focusing with exponent p > 1. In optics, the z variable in Eq. ([1]) is normalized 
by 2Ldiff, where is the diffraction (Rayleigh) length and the Xj variables are 

normalized by the input beam radius. 

We delineate several cases for the NLS JT]): 

4 

< p— 1 < -, the subcritical case, 
a 

4 

p — 1 = — , the critical case, 
a 

4 

p — 1 > -, the supercritical case. (2) 

In the subcritical case, the solitary waves A = e wz u v {yi) of the NLS (pQ) are stable, 
while in the critical and supercritical cases the solitary waves of the NLS ([1]) are 
unstable. The profile of a stable solitary wave experiences only minor changes under 
small perturbations as it propagates. On the other hand, unstable solitary waves can 
change dramatically due to the effect of an infinitesimal perturbation. For the NLS ([I]), 
unstable solitary waves either collapse after propagating a finite distance, or diffract as 
z goes to infinity [TJ [2] • 

Solitons have been thoroughly studied in view of their potential application in 
optical communications and switching devices (in nonlinear optics) or in quantum 
information science (in BEC). Recent advances in fabrication and experimental methods 
now make possible the realization of transparent materials with spatially varying, high 
contrast dielectric properties. Such materials have various all-optical signal processing 
applications in optical communications, see e.g. [31 0]. In this case, the solitons are 
usually called lattice solitons. Specifically, by a proper design of the dielectric properties 
of the medium, it may be possible to avoid the blowup/diffraction instability in the 
critical and supercritical cases and to obtain stable propagation of laser beams in 
those structures [SI El [3, [8] . Thus, there is considerable interest in understanding the 
propagation of light in modulated media. 
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Most studies of such media have considered linear lattices (potentials). In this case, 
the equation of propagation is 



where x; at = (xi, . . . , Xd lat ) are the lattice coordinates, 1 < d[ at < d is the lattice 
dimension and 1/N is the characteristic length-scale of change of the lattice. For 
example, if the lattice is periodic, then N is the lattice period. In the context of 
nonlinear optics, linear potentials are created by modulating the linear refractive index 
n in space. If the modulation/potential is periodic, such structures are called waveguide 
arrays or photonic lattices. In the context of BEC, the corresponding Gross-Pitaevskii 
equation accounts for the interaction of the atoms with a magnetic trap or, in the case of 
a periodic optical lattice, with interfering laser beams, see [HI [H)] and references therein. 

Solitary waves of the NLS ([3]) with a general linear potential were studied in [TT|[T2]. 
to name a few of the earlier studies. Recently, many studies considered periodic 
potentials. Theoretical and numerical studies of solitons of the NLS/GP equation were 
done for a periodic potential in one [131 O CESl CES] , two [13 [HI [19] and three [201 OH] 
dimensions. Experimental realization of these solitons were obtained in one-dimensional 
waveguide arrays [22] and in two-dimensional optically induced photonic lattices in 
photorefractive media [23l EH [25J, [26j [27]. Some studies also involved lattices whose 
dimensionality is smaller than the spatial dimension, i.e., d lat < d (see e.g. [8j [28]) and 
in media with a quintic nonlinearity (see [29] and references therein). 

Generally speaking, it was found that for some lattice types and propagation 
constants z/, the lattice can prevent the collapse and stabilize the solitons in the critical 
and supercritical cases. However, the possibility that these stable solitons can collapse 
under a sufficiently large perturbation was not mentioned in previous studies. 

A detailed study of stability (and collapse) of solitons in a nonlinear lattice, i.e., 



was done in j30[ [31]. In these studies it was shown that the soliton profile and 
(in) stability properties strongly depend on whether it is wider than, of the same order 
of, or narrower than the lattice period. Specifically, it has been shown that the same 
nonlinear lattice may stabilize beams of a certain width while destabilizing beams of 
a different width. Hence, any study of the stability of lattice solitons should take into 
account the (relative) soliton width. 

In this paper, we conduct a systematic study of the stability and instability 
dynamics of solitons in linear lattices which are narrow with respect to the lattice 
period. The fact that the solitons are narrow imply that there is a small non- 
dimensional parameter N, see Eq. ([6]). This allows us to employ perturbation methods 
and to compute the soliton profile and related quantities (soliton power, perturbed 
zero-eigenvalues Aq^\ see below) asymptotically. 

In nonlinear optics, typical lattice periods are of the order of several microns and 
typical input beam sizes are not smaller than this period [221 [2S1 [321 [331 E]- Hence, 
typically, the input beam sizes are not small compared with the lattice period. However, 




(3) 




(4) 
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if the beam undergoes collapse, the beam can become much narrower than the lattice 
period. In BEC, the standard magnetic traps are significantly wider than the size of the 
condensate. Hence, the narrow beams limit is of physical relevance. From a theoretical 
point of view, the limit of narrow beams corresponds to the semi-classical limit of the 
nonlinear Schrodinger equation 

ihA z (z, x) = -h 2 V 2 A - \A\*~ X A + V(x lat )A, h^O, (5) 

see e.g., [HI 154"] . Moreover, as discussed in Section El in many cases, the results for 
narrow beams hold also for beams of 0(1) width. 

The paper is organized as follows: In Section [2j we present various physical models 
in nonlinear optics and in BEC where Eq. ([3]) arises. In Section [31 the equation for lattice 
soliton is derived. It is shown that the soliton width is given by a single parameter 
N 

N= — <1, v = v + V(0), (6) 

where V(0) is the potential at the soliton center. Therefore, the limit v — > oo analyzed 
in [35], and the limit iV — > analyzed in [31], are in fact the same limit. It is well 
known that narrow solitons of a periodic lattice are found deep inside the "semi-infinite 
gap" of the linear problem, away from the first band of the allowed solutions [TB], i.e., 
for v — > oo. Indeed, in this case N — * 0. However, from this argument it is not clear 
how large should v be in order for the soliton to be narrow. This information is given 
by the parameter N, which is thus, a more informative parameter than the propagation 
constant v. Moreover, the parameter N includes also the effect of the lattice strength on 
the width and reflects the fact that as V(0) increases, the beam confinement increases, 
hence the beam becomes narrower 0. 

In Section [3j we also use perturbation analysis to calculate the profile of narrow 
lattice solitons for any dimension d, lattice dimensionality d\ at and nonlinearity exponent 
p. As can be expected, this calculation shows that the soliton profile depends only on 
the local properties of the lattice, rather than on the full lattice structure. Hence, our 
study is relevant to any slowly varying lattice, regardless of its long-scale properties. To 
simplify the notation, we mostly consider lattices that are aligned in the directions of 
the Cartesian axes. In this case, the lattice can be expanded as 

(dlat \ 
N2 J2v 2 j + ^(N 4 )J ■ (7) 

Our results are valid, however, to any linear lattice, see Remark 13.11 

In Section HI we analyze the stability of narrow lattice solitons. We first present 
the two conditions for stability of lattice solitons in Theorem 14.11 The first condition, 
known as the Vakhitov-Kolokolov condition [36] or the slope condition [37], is that 
the power (or L2 norm) of the soliton should increase with v. Using the results of the 
perturbation analysis, we show in Section H~Tl that to leading order, the power of a narrow 
lattice soliton is equivalent to the power of a soliton in a homogeneous medium, and 



% Note, however, that expression (128 j) for the beam relative width is only valid for narrow beams. 



that the change in the power due to the lattice scales as N 2 . [§| In particular, the lattice 
causes the power to decrease (increase) for lattice solitons centered at a lattice minimum 
(maximum). In addition, the power curve slope is more positive (negative) for lattice 
solitons centered at a lattice minimum (maximum). Since in a homogeneous medium 
the slope has an 0(1) magnitude in the subcritical and supercritical cases, the small 
change of the slope by the lattice does not affect the sign of the slope. Accordingly, the 
slope condition remains satisfied in the subcritical case and violated in the supercritical 
case. In the critical case, the slope in a homogeneous medium is zero. As a result, 
the 0(N 2 ) change in the power by the lattice leads to a positive (negative) slope for 
lattice solitons centered at a lattice minimum (maximum). Hence, the slope condition 
is satisfied for narrow lattice solitons centered at a lattice minimum, but is "even more" 
violated for lattice solitons centered at a lattice maximum. 

The second condition for stability of narrow lattice solitons is the spectral 
condition [39], and it involves the number of negative eigenvalues of the linearized 
operator L+l, see Eq. ( 134"1) . In Section l4~2l we first show that the spectral condition is 
violated if and only if the lattice causes some of the zero eigenvalues of the homogeneous 
medium linearized operator L+ )V (see Eq. (I42p ) to become negative. Then, we use a 
perturbation analysis to show that the values of the perturbed zero eigenvalues Xqj are 
given by 

{N) = f 5v J3 N 2 + 0(A> 4 ), j = d lat , 
° J 0, j = d lat + l,...,d, 



where 



. p(2 - d) + 2 + d 
o = — 



p — 1 

see Lemma 14.21 This calculation shows that the eigenvalues become positive (negative) 
for solitons centered at a lattice minimum (maximum). Hence, the spectral condition 
is satisfied (violated) for solitons centered at a lattice minimum (maximum). This 
calculation generalizes the result of Oh in the one-dimensional cubic case [TT] to any 
dimension d, any lattice dimension di at and any nonlinearity exponent p. 

In order to test the validity of the analytical formula for Aqj , we also compute these 
eigenvalues numerically. For d > 2, the matrix that represents the linearized operator 
L+l is very large. As a result, standard numerical schemes (e.g., Matlab's eig or eigs) 
usually fail to compute its eigenvalues. In order to overcome this numerical difficulty, 



we use a numerical scheme which is based on the Arnoldi algorithm, see Appendix C 



While in this study we "only" use this scheme to verify the validity of the analytical 
approximation of the eigenvalue, we note that in the case of non-narrow lattice solitons, 
the eigenvalue cannot be computed analytically, and the only way to check the spectral 
condition is numerically. Moreover, this numerical scheme can be used in similar 
eigenvalue problems in which large matrices are involved. 

§ For comparison, the change in the power due to a nonlinear lattice is 0{N 2 ) in the subcritical and 
supercritical cases but 0(N A ) in the critical case [5CT1I55] , 
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lattice minimum 


lattice maximum 


Subcritical 


Stability 


Instability^ (drift) 


Critical 


Stability 


Instability*'^ (blowup+drift) 


Supercritical 


Instability* (blowup) 


Instability*'"!" (blowup+drift) 



Table 1. Stability of narrow lattice solitons. Condition leading to instability is marked 
by * for a failure to satisfy the slope condition and by ' for a failure to satisfy the 
spectral condition. In the case of instability, its dynamics is indicated in parentheses. 

Combining the results of Sections 14.11 and 14.21 we show in Section 14.31 
(Proposition 14.21) that in the subcritical and critical cases, narrow lattice solitons are 
stable when centered at a lattice minimum, and unstable when centered at a lattice 
maximum or at a saddle point. In the supercritical case, narrow lattice solitons are 
unstable at both lattice maxima and minima. 

Proposition 14.21 specifies when the two conditions for stability are violated. It does 
not, however, describe the resulting instability dynamics. The relations between the 
condition which is violated and the instability dynamics were observed in [301 131] for a 
nonlinear lattice and in [40J for a linear delta-potential to be as follows: 

(i) if the slope is negative, the soliton width can undergo significant changes. In the 
critical and supercritical cases, this width instability can result in collapse. In 
the subcritical case, this width instability can "only" result in a "finite-width" 
instability, i.e., the soliton width can decrease substantially, but not to zero. 

(ii) When the spectral condition is violated, the solitons undergo a drift instability, 
i.e., the soliton drifts away from the lattice maximum towards the nearest lattice 
minimum. 

(iii) When both conditions for stability are violated, a combination of a width instability 
and a drift instability can be observed. 

In the case of narrow lattice solitons, the slope is always positive in the subcritical case. 
Hence, the instability due to a negative slope is a blowup instability and not a "finite- 
width" instability. Furthermore, in Section l4T4"l we prove that when the spectral condition 
is violated (i.e., if the soliton is centered at a lattice maximum or saddle point), narrow 
lattice solitons undergo a drift instability, i.e., they move away from their initial location 
at an exponential drift-rate. In contrast, solitons centered near a lattice minimum (for 
which the spectral condition is satisfied) undergo small oscillations around the lattice 
minimum. The above observations on the condition leading to instability and the type 
of instability dynamics are summarized in Table [U 

In Section [5J we study the dynamics of solitons in the two cases where the small 
effect of the lattice changes the stability. As observed in [30, 31J, in such cases, it 
is important to study both stability and instability quantitatively. In Section 15.11 we 
discuss the strength of the stabilization induced by the lattice for solitons centered at a 
lattice minimum in the critical case. To do so, we use the concept of the stability region, 
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i.e., the region in function space of initial conditions around the soliton profile that lead 
to a stable propagation. As in the case of a nonlinear lattice [30|, [3Tj. our results indicate 
that the 0(N 2 ) small slope of the power curve implies that the stability region is 0(N 2 ) 
small [jj. Therefore, although the two conditions for stability are satisfied, these solitons 
can become unstable under extremely small perturbations. Practically, this means that 
in the critical case, "mathematically" stable solutions can be "physically" unstable i.e., 
become unstable under typical physical perturbations. We illustrate these results using 
two standard types of lattices: A sinusoidal potential, which is typical in photorefractive 
materials [23l [25] and in BEC [4Tj and a Kronig-Penney step lattice (periodic array of 
finite potential wells) [42J, which is typical for manufactured slab waveguide arrays, 
see e.g., [T3| [22| [33] . We study numerically the stability of solitons under random 
perturbations that either increase or decrease the total power of the soliton and observe 
that narrow lattice solitons are "mathematically" stable but "physically" unstable. The 
stability is particularly weak for Kronig-Penney lattice solitons, for which the slope is 
exponentially small. In addition, we observe that when the perturbation is sufficiently 
"non-small" , both the sinusoidal and KP (stable) lattice solitons can undergo a blowup 
instability. This shows that in the absence of translation invariance, stability and blowup 
can coexist in NLS equations [30 | I3TI 143] . 

In Section 15. 2\ we show that the opposite scenario is also possible, i.e., 
"mathematically unstable" solitons can be "physically stable". This occurs for 
subcritical narrow lattice solitons centered at a lattice maximum, which are unstable 
due to a violation of the spectral condition (Proposition 14.21) . We show that the drift 

1/2 

rate is exponential in y — Aq^J • Therefore, narrow solitons, for which Aq^ is 0(N 2 ) 
small, experience very slow drift and can thus be "stable" for the distances/times in 
experimental setups. In particular, we observe that the Kronig-Penney lattice soliton 
drifts much more slowly than the sinusoidal lattice soliton of the same width. Section [6] 
concludes with some concluding remarks. 



2. Physical models 

We consider the d dimensional NLS equation (j3J) with a linear lattice in di at dimensions 
(1 < diat < d). This model describes numerous physical configurations. For example, 
beam propagation in a Kerr slab waveguide with a lattice is modeled by 

iA z (z,x) = -A xx - \A\ 2 A + V(Nx)A. (8) 

In this case, p = 3, d = d\ at = 1, x = xi at = x and V = V(Nx), see e.g., [HJ, [151 122] . 
Beam propagation in bulk Kerr medium with a two-dimensional lattice is modeled by 

iA z (z,x,y) = -V 2 A- \A\ 2 A + VA. (9) 

|| In the case of a nonlinear lattice, the slope, hence the size of the stability region, is 0(N 4 ) small, 
implying an even weaker stability [301 138j . 
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In this case, p = 3, d = 2 and x = (x,y). If V — V(Nx, Ny), then di at = 2, and 
*iat = (x,y), see e.g., pH [191 E]; if V" = V^iVx) then d iat = 1, and x Zat = x. In the 
latter case, the dimension of the lattice d\ at is smaller by one from the dimension of the 
transverse space d, see e.g., [51 l3T| [28]. 

Propagation of ultrashort pulses in a slab waveguide is modeled by 

iA z {z, x, t) = -A xx + faA tt - \A\ 2 A + V{Nx)A, (10) 

where fa is the group velocity dispersion (GVD) parameter. In the case of anomalous 
dispersion (fa < 0), the time coordinate t is effectively an additional transverse 
dimension. Then, Eq. (flQj) corresponds to Eq. (J3]) with p = 3, d = 2, x = (x,t), 
di a t = 1 and x; at = x, so the dimension of the lattice di a t is smaller by one from the 
dimension of the transverse space d, see e.g., [5J,[3T]. Similarly, propagation of ultrashort 
pulses in a 2D optical lattice is modeled by 

iA z {z, x, y, t) = -V 2 A + (3A tt - \A\ 2 A + VA, (11) 

which for (3 < corresponds to Eq. (j5j) with p — 3, d = 3 and x = (x,y,t). If 
V = V(Nx,Ny), then d lat = 2, and x /at = (x,y) [2TJ; if V = V(Nx) then <i/ at = 1, and 

*lat = X. 

The linear lattice V in Eq. ([3]) varies in the transverse coordinates but not in z. In 
some applications, the lattice varies in the direction of propagation z. Such problems, 
however, will not be studied in this paper. 

Eq. ([3]) also models the dynamics of Bose-Einstein condensates (BEC) with a 
negative scattering length. In this case, z is replaced with t. In BEC, typically 
x = (x,y,z), i.e., d = 3, but under certain conditions the cases d — 1 and d = 2 
are also of physical interest, see e.g., [151 SB]- The exponent p is usually equal to 3 
but can also be equal to 5, see [29] and references therein. In the BEC context, both a 
parabolic potential and a periodic potential appear in the experimental setups [9J. 

3. Narrow lattice solitons 

We look for lattice solitons, which are of solutions of Eq. ([3|) of the form 

A(z,x) = e ivz u^(x), u>0, (12) 

where Uu^ is the solution of 

V 2 «f)(x) + (u^f -[u+ V(Nx lat )] U W = o. (13) 

We consider lattices which are symmetric with respect to a critical point x|°] of the 
lattice V. Hence, the soliton maximal amplitude is attained at x z °] [47] . The boundary 
conditions for Eq. ( fl3|) are Vu„ (x|°|) = and ui N \oo) = 0. Without loss of generality, 

We study solutions of Eq. (fl3l) which are narrow with respect to the lattice 
characteristic length-scale. A priori, the relative width of a lattice depends on the 
lattice strength, the lattice period (or characteristic length) 1/N and the propagation 
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constant v. We now show that in the case of narrow solitons, one can rescale Eq. (1131) 
to a form where the relative width of the beam is given by a single parameter N. In 
order to achieve that, we define 

V = v + V{0), N = N/^j, «W(x)=^(V^). (14) 

Then, Eq. ffTB"]) becomes 

V 2 M*)+4- [l + V{N5i lat ))u R = 0, Vujv(O), u#(oo) = 0, (15) 

where 

- r ~ r f/(Kr ~ s Vm lat ) - V(0) 

X = Xiat = y/rjxiat, V (NyL lat ) = . (16) 

When N <C 1, we can expand the solution of Eq. (115]) as a power series of N 2 , i.e., 

u#(x) = W(|x|) + N 2 g(±) + 0(N 4 ), (17) 

where U is the positive, radially-symmetric ground-state solution of 

V 2 W(|x|) +U P -U = 0. (18) 

Similarly, since V"(0) = and VV"(0) = 0, the potential V(N~ki at ) can be expanded for 
N < 1 as 

^(iVx Zat ) = iV 2 V 2 (x^) + G(iV 4 ), (19) 

where 

V2(Xlat) = }^ V jk XjX k , V jk = - (20) 



2 dy 3 dy k 



Ylat=0 



is the first non-vanishing term in the Taylor expansion of V which represents the local 
curvature of the lattice at the soliton center. In particular, V^x^t) > (< 0) for lattice 
solitons centered at a lattice minimum (maximum). 

Remark 3.1 In order to simplify the presentation, we assume that the principle axes 
of the lattice identify with the Cartesian axes {e 1: . . . , e^ at }. In this case, Vj k = for 
./ :' /■•• 

dlat 

3=1 

and 

(dlat \ 
N2 J2v 2 3+°(N 4 )j , (22) 

see Eq. AC.3\) . However, all our results can be immediately generalized to the case when 
the lattice is not aligned along the cartesian axes as follows. Since Vj k = v k j, there exists 
a basis of vectors {ex, . . . , £d lat } such that if 5q ai = Ylj=l a j^j then 

V 2 (± lat ) = Y,^a 2 . (23) 

3=1 
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Therefore, in order to apply our results to the lattice Kffl^l , one needs to replace Xj by a>j 
andvjj byujj. See e.g., Reraark WT^ and Remark \^A , 



Using a perturbation analysis similar to the one used in [301 EE], we show that 
Lemma 3.1 The solution of Eq. (T5\) for N <C 1 is given by 

u*(x) = U(\5t\) - N'L- 1 (v 2 {% at )U) + 0(N A ), (24) 

where U is given by Eq. [T8\) . V 2 is given by Eq. WIS) and 

L + = - V| - pU^ 1 + 1. (25) 



Proof: See Appendix A 



u 



In the original variables, the expansion fl24j) becomes 
> % x ) = {u + V(0))p^ \u{y/v + V(0)\x\) + 0(N 2 )} = W„(|x|) + 0(N 2 ), (26) 



where = rjp- 1 U( y /rj\x.\) is the solution of 

VX(N)+Z<-^ = 0. (27) 
This expansion shows that: 

(i) To leading order, a (rescaled) narrow lattice soliton is given by the rescaled 
homogeneous medium soliton U. 

(ii) The deviation of a narrow lattice soliton from U is 0(N 2 ) small, even if the 
lattice has 0(1) variations. 

The above results also show that the soliton relative width is given by a single 
parameter A^: 

Proposition 3.1 Lattice solitons are narrow with respect to the lattice period if 
N N 

N = — = - < 1. 28 

Vv \/v + V(0) 

In this case, 

^ soliton width 
lattice period 

Proof: By Eq. (f2"4"|l . when JVC 1, then Ufj has 0(1) width in x. Hence, by Eq. f T26|) . 

the width of ui N \x) in x is 0(l/+/r}). Since that the lattice length-scale/period is 1/N, 

then the relative width of the soliton is given by N. □ 

We emphasize that the expansion (124"]) applies to all types of lattices so long as 

v jjN 2 «C 1. Specifically, for a strong periodic lattice (V ^> 1), for which the linear 

coupling between adjacent lattice sites is weak, the result (1261) is still valid provided 

_i 

that A^ is small enough, i.e., for A" -C v^ 2 . In that case, the solution (j2BJ) is the 
continuous analog of the discrete solitons of the DNLS model 
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3.1. Effect of lattice type 

Lemma 13.11 shows that the effect of the lattice depends on whether V 2 ^ or V 2 = 0. 
When V 2 ^ 0, then the lattice effect is 0(N 2 ). This case corresponds to a parabolic 
lattice, a sinusoidal lattice etc. However, when V 2 = 0, then g(x) = and the next-order 
term in the expansion (fTTl) must be considered. In particular, in the special case of a 
Kronig- Penney step lattice (see, e.g., Eq. fl38|) ). all derivatives of V at the soliton center 
x|°] = vanish. Therefore, the difference between and U will be exponentially small. 

3.2. Effect of lattice inhomogeneity on soliton profile 

In order to calculate the effect of the lattice on the soliton profile, we note that 

L- l (x)U) =4S(|x|) + Q(|x|), (30) 
where S and Q are radial functions which are the solutions of 

L+S - ts' = U, L+Q = 2S. (31) 

r 

Indeed, applying the operator L + to the right-hand-side of Eq. (1301 gives 
L + (4S-(|x|) + Q(|x|)) = (-Vl + l-pW- 1 ) (x*S(r) + Q(r)) 

= x) U + S{r) - ts'{r)\ - pS(r)+L+Q(r)) =xft. 



=0 



=u 

t2\ 



Therefore, the 0(N ) correction to the soliton profile due to the lattice (T2T1) is given by, 
see Eq. ([21]), 

d-lat dlat 



-U~ -N'L- 1 (V 2 (5c lat )U) = -iV 2 (5(|x|)^ v ]3 ~x) +g(|x|)£; % ). (32) 



anisotropic isotropic 

Thus, the variation of the lattice in the direction Xj has an isotropic effect through 
<5(|x|) and an anisotropic effect in the direction Xj through 5(1x1). 

Remark 3.2 IfV 2 is given by the lattice ( TJS| ), then, the 0(N 2 ) correction to the soliton 
profile is given by 

(dlat d iat \ 

5(1x1)^^ + ^(1x1)^^, . (33) 
3=1 j=l J 

4. Stability and instability of lattice solitons 

Eq. fl26|) implies that narrow lattice solitons ui are positive. The conditions 
for stability and instability of positive lattice solitons are as follows ([50]. and see 
also [191 [371 [39J ) : 
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Theorem 4.1 Letu^ be a positive solution of Eq. (T7^). letvj) N ^ = / ^ui^J ^ x ^ e ^ e 

power of ui N \ and let n_(412) fre i/ie number of negative eigenvalues of the linearized 
operator 

45 = -V 2 + i/ - p(«^(x)) p - 1 + K(iVx^). (34) 
Then, the lattice soliton A(z,x) = e wz Uv^ (x) zs 

(%) an orbitally stable solution of Eq. |5|) if 

(a) d v vi N ^ > (slope condition), and 

(b) n~{L+l) = 1 (spectral condition). 

(ii) an orbitally unstable solution of Eq. (Ej) 

(aj d v Vu N ^ < 0, or 
(b) n_(L<$) > 1. 

In what follows, we use the expansion (124]) to determine whether the two conditions in 
Theorem 14.11 are satisfied ,and consequently determine the stability of narrow lattice 
solitons. 

4-1. Slope condition 

We can use the expansion ( 1241) to calculate the power of narrow lattice solitons: 
Lemma 4.1 The power of narrow lattice solitons (N <C 1) is given by 

(dint \ 
V u=1 - C V N 2 J2 «ii + °(^) J , (35) 

where V v =\ = J \U\ 2 d5t, U is the positive solution of Eq. / figj) . and 

Cy = 2P ~l + dp ~ d ! \x\Wdx (36) 



2d(p - 1 

is a constant independent of N and v. 



Proof: See Appendix B 



Eq. (1351) shows that, in a similar manner to its effect on the soliton profile, when 
V2 ^ (e.g., in the case of a sinusoidal or parabolic lattice), the lattice has an 0(N 2 ) 
small effect on the soliton power, even if the lattice itself is not weak. In light of 
Section 13.11 in the case of a Kronig- Penney step lattice, the effect of the lattice on the 
power is exponentially small in N. 

^From Eq. (I3"6"j) it also follows that Cy > for p > 1 + ^r^. In particular, for p = 3, 
C v = \j r 2 U 2 d± > 0. Thus, 

Corollary 4.1 If V2 ^ ; the lattice causes the power to decrease (increase) for lattice 
solitons centered at a lattice minimum (maximum) for anyp > l+d^; an d in particular, 
for a Kerr nonlinearity p = 3. 
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In order to demonstrate the result of Lemma 14.11 we solve Eq. (1151) numerically 
with d = di a t = 2, x = x iat = (x, y) and p = 3. For convenience, the numerical results 
shown here are presented for rj = 1 (so that N = N, Uv = u^) and for V(0) = (so 
that V — V). We study two different two-dimensional lattices with a periodic square 
topology: 

(i) A 2D sinusoidal lattice given by 

V(Nx, Ny) = ±- (sin 2 (vriVx) + sm 2 {irNy)) . (37) 

(ii) A 2D Kronig-Penney lattice that consists of an array of primitive cells of size 
[-1/N l/N] x [-1/N l/N], each consisting of circular waveguide with abrupt 
index change between and 1, i.e., 

I ±1, otherwise. 

In both cases, the plus/minus sign corresponds to a lattice with a minimum/maximum 
at x c = 0, respectively. The parameters of these lattices were chosen so that both 
lattices have a period 1, mean value 1/2 and vary from to ±1. The lattices are shown 
in Fig. [T] for a lattice with a minimum at x c = 0. Note that both lattices are anisotropic 
in r = \J x 2 + y 2 , and thus, require a full 2-dimensional treatment. Moreover, since the 
2D cubic NLS is critical, V v= \ = V cr — 11-7, where V cr is the critical power for collapse 
in a homogeneous Kerr medium. 

In Fig. [21 we show the power of narrow lattice solitons centered at a lattice 
minimum for both lattices. For < N < 0.1 there is good agreement between the 
numerically calculated value of the power of the sinusoidal lattice solitons and the 
analytical approximation 

diat 

p( N ) = Vv=x _ C yN 2 v n + £ 11.7 - 6.94 • 2(2n 2 )N £ 11.7 - 273.8A> 2 , (39) 

i=i 

which is derived from Lemma ETTl In particular, the effect of the lattice on the power of 
the narrow lattice solitons is much more pronounced in the case of a sinusoidal lattice 
than in the case of a Kronig-Penney lattice. 



The sign of the slope follows directly from Eq. (135]) : 

Corollary 4.2 Let N <C 1. Then, the slope dyPl^ is positive in the subcritical case 
(p < 1 + 4:/d) and negative in the supercritical case (p > 1 + A/d). In the critical 
case (p = 1 + A/d), the slope is positive for narrow lattice solitons centered at a lattice 
minimum and negative for narrow lattice solitons centered at a lattice maximum. 

% The agreement between the analytic result ([55)1 and the numerics is good "only" for relatively small 
values of N because of the large curvature (X)j=i v jj = 4t 2 ) of the lattice which translates into a large 
coefficient of the N 2 term in Eq. (j39|) . Indeed, we verified that for smaller values of Ylj=i v jji ^ ne 
agreement between the analytic result (|35|) and the numerics extends to larger values of N. 
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Figure 1. (a) sinusoidal lattice ([57)1 (b) Kronig-Penney lattice ([55]) . 




N 



Figure 2. Relative power deviation from V v= \ — V C r for narrow sinusoidal (dots) and 
Kronig-Penney (dashes) lattice solitons centered at a lattice minimum. The analytical 
prediction ([39| for the sinusoidal lattice is shown by a solid line. 



Proof: In the subcritical and supercritical cases, the slope is given by 

N 2 



d„vl N) 



4-d(p-l) 



v v=l + o 



4-d(p-l) 

-1) 



o 



is + V(0)J 

N 2 



(40) 



4-d(p-l) -i 



v + V(0) 

Therefore, in the subcritical case, the slope is positive while in the supercritical case, 
the slope is negative. Note that in these cases, the lattice does not affect the sign of the 
slope. 

In the critical case, the first term in Eq. (HOj) vanishes and the slope is determined 
by the 0(N 2 ) correction in Eq. (I35I) . i.e., 

dint 



QA7-2 a '«i 



2CY 



N 



(41) 
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where we also used Eq. (JHJ). By Eq. (}3T)j) . in the critical case Cy — h I f 2 U 2 d~k > 0, 
which completes the proof. □. 

We thus conclude that although the lattice has a small effect on the profile of narrow 
lattice solitons, in the critical case, this small effect determines the sign of the power 
slope and hence, the stability (but see Section l5J\) . 

4-2. Spectral condition 

As noted is Section HJ lattice solitons are stable only if in addition to the slope condition, 
they also satisfy the spectral condition. In the absence of a lattice (i.e., for V = 0), the 
linearized operator L+l reduces to L +)V which is given by 

L +>l/ = -V 2 -pUr 1 + ^, (42) 

where U v = u~U{y/u\y.\) and U is given by Eq. ( fT8l) . The spectrum of L +j „ consists 
of 



(i) A negative eigenvalue X min and a corresponding even and positive eigenfunction 
fu,min- In [EE], Oh shows that for d — 1 and p = 3, X min = —3u and f v ,m%n — ■ 
More generally, we observe that for any value of p and d, 

1 s+l 

X m in = - J (P ~ 1) (P + 3) U, f v , m in = U v 2 . 

(ii) A zero eigenvalue Ao of multiplicity d with the corresponding eigenf unctions 

f»Ax) = ^ = rrK(\x\), j = i,...,d. (43) 

OXj |x| 

(iii) A positive continuous spectrum [u,oo). 

Thus, in a homogeneous medium the spectral condition is satisfied. In the presence 
of a linear lattice, the perturbed smallest eigenvalue A^ remains negative. The 
continuous spectrum develops a band structure, but remains positive. Moreover, for 
diat < j < d, the jth perturbed zero eigenvalue remains at zero with the corresponding 

r) ( N ) (N) 

eigenfunction ^ . Therefore, L + J, can attain more than one negative eigenvalue only 
if at least one A^ becomes negative for 1 < j < d\ at |30j. Thus, in order to check if the 
spectral condition is satisfied, we only need to compute the sign of A^ for 1 < j < d iat . 

For d — 1, p — 3 and a slowly varying parabolic potential, the value of the perturbed 
zero eigenvalue A^ = A^ was computed by Oh [IT] : El 

Aq^ = 3vjjN 2 + 0(N 3 ). (44) 

A more general result on the value and sign of Aq^ in the presence of a linear lattice for 
d > 2 is not known to us. We now give an asymptotic formula for XqO for narrow lattice 
solitons which generalizes of the result of Oh to any dimension d, lattice dimension d\ a t 
and nonlinearity p: 

+ The formula given in contains a minor error, since in pp. 29 of the Li norm of U was used 
instead of the Li norm of U' . 
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Lemma 4.2 Let V be given by Eq. (dJ|) ; or equivalently, let V2 be given by Eq. 

N 



and let N 1 . Then, the perturbed zero eigenvalues of the operator L+J are given 



by 



where 



(JV) _ 1 5 Vjj N 2 + (9(AT 4 ), j = 1, . . . , dfct, 

0J ~ i 0, j = d lat + i,...,d, (45j 



tf = P (2-d) + 2 + d 
p — 1 



Proof: See |Appendix C 



Remark 4.1 IfV has the general form j[W\) . then, Eq. fijlfy becomes 

^ = ^N 2 S + 0(N A ), j = l,..., d lat , (47) 



and Eq. (J/.6) remains unchanged. 
Proposition 4.1 Let 

1<P, d=l,2 

1 < p < g|, d> 2 ■ l48J 

Then, the spectral condition is satisfied for narrow lattice solitons centered at a lattice 
minimum, and violated for narrow lattice solitons centered at a lattice maximum. 



Proof: It is easy to verify that 8 > if and only if p satisfies condition (1481) . Thus, 
Lemma [4.21 shows that 

s gn( X aj) = sgniyjj). 

Consequently, the operator L^} has one negative eigenvalue (A^ > 0) for a narrow 
lattice soliton centered at a lattice minimum (vjj > 0) and more than one negative 
eigenvalue (A^ < 0) for a narrow lattice soliton centered at a lattice maximum 
(% < 0). □• 



We note that values of p for which condition (1481) is satisfied include all the 
physically relevant cases of d = 1, 2, 3 and p = 3, 5. 

To demonstrate the results of Lemma 14.21 we consider the case of d = d\ a t = 2, 
p = 3 and the lattice ([37}. By Eq. (SSD , 

Aff = Ag } = 2v 3J N 2 = ±(2n) 2 N 2 . (49) 

In order to confirm the validity of the expansion (l4"9l) . we compute the eigenvalues of 
the discretized operator L^J for the lattice (l37j) . In general, for d > 2, computation of 
the eigenvalues of the discretized operator L+l (using, e.g., Matlab's eig or eigs) 
fails to give reliable solutions due to computer memory limitation. In order to 
overcome this limitation, we used an improved numerical scheme based on the Arnoldi 



algorithm (see Appendix D ). In Fig. [3] we see that indeed for IV < 1, the asymptotic 



expression (1491) for the eigenvalue is in agreement with its numerically calculated value. 
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N 



Figure 3. Eigenvalue Aq - (j = 1,2) of the operator j as a function of jV 
for the lattice (|37|) and a for soliton centered at a lattice minimum (left) and at 
a lattice maximum (right). For N <C 1, there is a good agreement between the 
numerically calculated eigenvalue of the discretized operator (dots) and the 

analytical approximation (|49|) (solid line). 



4-3. Stability results 

Now that we have determined when the slope and spectral conditions are satisfied, we 
can characterize the stability of narrow lattice solitons: 

Proposition 4.2 Let N <C 1, let ui N ^ be the solution of Eq. / f73j) . let p satisfy 
conditions and let V be given by Eq. Then, 

(i) If Uv is centered at a lattice maximum, then ui e il/z is unstable. 

(ii) If ui N ^ is centered at a lattice minimum, then u^\ tuz is stable in the subcritical 
and critical cases p < 1 + A/d, and unstable in the supercritical case p > 1 + 4/d. 

Proof: Instability of narrow lattice solitons centered at a lattice maximum follows 
from a violation of the spectral condition (Proposition 14. Xf) . For narrow lattice solitons 
centered at a lattice maximum the spectral condition is satisfied (Proposition 14. lj) and 
stability is determined by the slope condition. Hence, the stability in the subcritical 
and critical cases and instability in the supercritical case follow from Corollary 14.21 □. 

Proposition 14.21 refers only to solitons centered at a lattice minimum or maximum. 
In some cases (e.g., in studies of lattices with defects or surface/corner solitons |50j). 
lattice solitons can be centered at critical points of the lattice that are saddle points. 
In these cases, by Lemma 14.21 the narrow lattice solitons are unstable since the spectral 
condition is violated. 

4-4- Instability dynamics 

Proposition 14.21 specifies the conditions for which narrow lattice solitons are unstable. It 
does not, however, describe the instability dynamics that occur when those conditions 
are not met. As noted in the Introduction, in previous studies [301 EH HO] it was 
observed that if the slope is negative, the solitons undergo a width instability and when 
the spectral condition is violated, the solitons undergo a drift instability. 
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In the case of narrow lattice solitons we can prove that violation of the spectral 
condition results in a drift instability by monitoring the dynamics of the soliton center 
of mass: 

Lemma 4.3 Let (xj) be the center of mass in the Xj coordinate, i.e., 

f xh I -^4 1 d'x. , 
= rliL • (5°) 



Then, 




(xj{0)} cos(fiz) + sm(Os), Vjj > 

(xj(0)) cosh(^) + sinh(n^), Vjj < 



where 



Proof: See Appendix E 



(51) 



n = 2N^dri\ Vjj \, (52) 
and Vjj defined in Eq. K21\). 



Thus, if Vjj > 0, the center of mass (xj) oscillates around the lattice minimum. On 
the other hand, if Vjj < 0, the center of mass moves away from the lattice maximum at 
an exponential rate. This shows, in particular, that a soliton centered at a saddle point 
is stable in the directions in which it is centered at a lattice minimum and undergoes a 
drift instability in the directions in which it is centered at a lattice maximum. 



5. Quantitative study of stability 

As noted, the lattice has a small 0(N 2 ) effect on the slope and on the value of the 
perturbed near zero-eigenvalues of Nevertheless, this small effect changes the 

stability of solitons centered at a lattice maximum (which became unstable) and of 
solitons centered at a lattice minimum in the critical case (which become stable). As 
pointed out in [301 EI] , when a small effect changes the stability, stability and instability 
needs also to be studied quantitatively. 



5.1. "Mathematical" stability vs. "physical" stability 

Let us first consider narrow lattice solitons centered at a lattice minimum in the critical 
case. In this case, according to Proposition 14.21 the solitons are stable. However, as 
was shown in [301 EI], satisfying the "mathematical" conditions for stability does not 
necessarily "prevent" the development of instabilities due to small perturbations. In 
order to understand how this can happen, we recall that Theorem 14. II ensures that there 
is a stability region in the function space of initial conditions around the soliton profile for 
which the solution remains stable. However, it does not say how large this stability region 
is. If the stability region is very narrow, the solution is only stable under extremely small 
perturbations. In this case, it is "mathematically" stable but "physically unstable", 
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i.e., it can become unstable under perturbations present in an experimental setup. If, 
on the other hand, it is also stable under perturbations comparable in magnitude to 
perturbations in actual physical setups, one can say that it is also "physically stable" . 

The distinction between "mathematical stability" and "physical stability" is only 
important in the critical case where, in the absence of the lattice, the slope is zero. 
Then, the slope (VK) condition shows that these solitons are unstable and indeed, an 
arbitrarily small perturbation can cause them either to undergo diffraction or to collapse. 
The effect of a linear lattice on narrow lattice solitons centered at a lattice minimum 
is to induce an 0(N 2 ) positive correction to the power slope which causes the slope 
(VK) condition to be satisfied and the solitons to become stable. As demonstrated for 
the first time in [301 EI] , the size of the stability region depends on the magnitude of 
the slope. This means that the transition between instability and stability is gradual 
rather than sharp, in the sense that as the soliton width N increases from zero, the 
magnitude of the slope grows from zero, hence the width of the stability region grows 
from zero. For example, in the case of a Kronig- Penney lattice, the power slope of narrow 
lattice solitons is exponentially small (see Section H~T1) . hence the stability region is also 
exponentially small. Therefore, narrow Kronig-Penney solitons are "mathematically" 
stable but "physically" unstable. On the other hand, in the case of a sinusoidal lattice, 
the stability region of the solitons is bigger, so that the sinusoidal lattice solitons can 
be also "physically" stable. 

In order to motivate the claims stated above, we first note that by definition (1251) 
of N, the slope d v Vv is proportional to df/Po . Thus, the slope with respect to the 
soliton width N can be viewed as a measure for the slope with respect to the propagation 
constant v. Second, we recall that the soliton profile Ufq is an attractor for NLS solutions. 
Therefore, small perturbations of the initial profile essentially lead to small oscillations 
of the soliton width along the propagation (see below). Thus, heuristically, we can view 
these width oscillations as a movement along the curve Vv- Such movement along 
the curve P„ was demonstrated e.g. in Fig. 6 of [ID]. Since the power is conserved, 
a large slope only allows for small changes of the soliton width (i.e., stability) while 
a small slope allows for larger changes of the soliton width and larger deviations from 
the initial state (i.e., instability). More generally, these arguments show that while the 
sign of the slope determines whether the solution is stable or not, the magnitude of 
the slope \d u Vi N \ corresponds to the size of the stability region. Hence, if the slope 
d v vl N ^ is positive but small, the stability induced by the lattice is weak. Therefore, if 
the perturbation applied to the narrow lattice soliton is large enough, the perturbation 
can "overcome" the stabilization and the solution will become unstable. 

A schematic illustration of the stability region in the critical function 
of the beam power V and the relative width N is shown in Fig. HI The stability 
region is centered around the lattice soliton power Vi N) = V cr - C V N 2 , see Eq. ([35]). 
By Eq. (j41j) and the above arguments, the size of the stability region depends on the 
propagation constant u, the period N and the lattice V(x) only through the parameter 
N, and is 0(N 2 ) small. Initial conditions to the left of the stability region undergo 
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a diffraction instability whereas initial conditions to the right of the stability region 
undergo a blowup instability. The separatrix between the stability region and blowup 
region can be estimated by the critical power for collapse in homogeneous medium V cr - 
Indeed, while the minimal power needed for collapse depends on the beam profile, for 
single-hump profiles such as ui N \ the minimal power needed for collapse is only slightly 
above V cr [SI]- 
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(b) 
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Figure 4. A schematic illustration of stability (shaded), diffraction instability and 
blowup instability regions as a function of the input beam width N and power V for 
narrow lattice solitons centered at a minimum of (a) a sinusoidal lattice and (b) a 
Kronig- Penney lattice. Dashed curve is vi N \ 



To illustrate these ideas numerically, we solve Eq. ([3]) for d = di at = 2 and p = 3, 
which correspond to the physical case of a 2D Kerr medium and N = 0.1 (i.e., narrow 
lattice solitons). Since this is the critical case, the lattice should have a dominant 
effect on the stability (see Proposition 14. 2\i . In order to demonstrate the difference 
between the stabilization by the sinusoidal lattice (1371) and by the Kronig-Penney 
lattice ( 1381) . we perform a series of numerical simulations with the initial condition 
A (x 1 y) = (1 + e • h(x,y))ui N \ Here is = r/ = I and h(x,y) is a random function 
which is uniformly distributed in [0,1] x [0,1]. Hence, the perturbation increases the 
power of the initial condition by the factor of ~ (1 +e) with respect to the power of the 
soliton Uv . We consider narrow solitons centered at a lattice minimum, hence they are 
"mathematically" stable, see Table [TJ 

We first note that in all the simulations in this Section, the center of mass of the 
beam, which is initially perturbed from the lattice minimum due to the random noise, 
remains small and close to the lattice minimum, in accordance with Lemma 14.31 

In Fig. [5](a), we show the solution for the Kronig-Penney lattice for various values 
of e > (i.e., when the noise increases the beam power) for < z < 70, i.e., over 
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Figure 5. Maximum intensity vs. propagation distance of narrow lattice solitons 
(N = 0.1) with power-increasing random perturbations for (a) Kronig-Penney 
lattice and (b) sinusoidal lattice (|3"T)) . Comparison of the dynamics for a sinusoidal 
lattice (solid) and Kronig-Penney lattices (dashed) is shown in (c) for e = 0.003. 



140 diffraction lengths. For e = 0.001 and 0.002, the solution undergoes a focusing- 
defocusing oscillations. When the initial perturbation is further increased (e = 0.003), 
the beam undergoes collapse. The abrupt change in the dynamics between e = 0.002 
and e = 0.003 can be understood by looking at the power of the beams. For the specific 
noise realizations in our simulations, the power of the initial condition was slightly 
below the critical power V cr for e = 0.001 and 0.002 and slightly above V cr for e = 0.003. 
Therefore, the beam undergoes collapse in the latter case. 

While an e = 0.003 perturbation to a Kronig-Penney lattice soliton leads to collapse, 
the same perturbation applied to a narrow sinusoidal lattice soliton only leads to small 
amplitude oscillations, see Fig. 0(b). When the perturbation is increased to e = 0.02 
the oscillations become stronger yet the solution does not collapse. Only when the 
perturbation is further increased to e = 0.035 the beam collapses in a finite distance. 
As in Fig. [3(a), we confirmed that for e = 0.003 and e = 0.02 the beam power is below 
P cr , while for e = 0.035 it is above P cr . 

These simulations confirm that although both lattice solitons are "mathematically" 
stable, sufficiently large perturbations can still cause these stable solitons to undergo 
collapse 0. This demonstrates that collapse and stability can co-exist, see also [i3| 138]. 
Moreover, these simulations also support the heuristic argument presented in Section [5711 
that the upper boundary of the stability region can be estimated by the critical power 
for collapse in homogeneous medium P cr . 

In Fig. EJ we show the solutions for e = —0.001 and e = —0.003 (i.e., when the 
noise decreases the beam power). The comparison between the two lattices for the same 
value of e shows that the stabilization by the sinusoidal lattice is much stronger than 
by a Kronig-Penney lattice. Additional simulations (data not shown) show that the 
difference between the stabilization by the two lattices becomes more pronounced as 
TV becomes smaller. Indeed, for a Kronig-Penney lattice, the boundaries of the lattice 
are located far in the soliton tail region. Thus, their presence can prevent broadening 

* Note that the typical perturbations in experimental setups are at least of few percents. 
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Figure 6. Maximum intensity vs. propagation distance of narrow (N = 0.1) lattice 
solitons with power-decreasing random perturbations for sinusoidal (solid) and Kronig- 
Penney lattices (dashed) with (a) e = —0.001 and (b) e = —0.003. 



only once the narrow beam has undergone significant broadening. On the other hand, 
a sinusoidal lattice acts at any position in the central region of the soliton, hence, it has 
a much more pronounced effect. 

The results shown in Figs. [5] and [6] confirm that Kronig- Penney lattice solitons are 
"physically unstable" (i.e., an extremely small stability region) whereas sinusoidal lattice 
solitons can be "physically stable" (not-so-small stability region). Indeed, a comparison 
between these two lattices for the same value of e shows that for narrow lattice solitons, 
the same perturbation leads to collapse in the case of a Kronig-Penney lattice but only 
to small oscillations and stable behaviour in the case of a sinusoidal lattice, see Fig. E{c) 
and Fig. [61 



5.2. "Mathematical" vs. "physical" instability 

We now consider narrow lattice solitons centered at a lattice maximum. According 
to Proposition 14.21 these solitons are unstable as they violate the spectral condition. 
Indeed, we showed that these solitons undergo a drift instability away from the lattice 
maximum. Since there is no drift for Aoj = 0, by continuity, the drift rate should be 
"small" for small negative values of Aq^. Indeed, combining Eqs. (1451) and (l5Tj) . one 
sees that for Vjj < 0, 



(xj(z)) ~ (^-(0)) cosh(fiz) + ^251 S inh(^), n = 2\j H^lAA. (53) 

Thus, if — Xq N ^ is small, the instability develops very slowly. In this case, the solitons 
are "mathematically" unstable but can be "physically stable", i.e., the instability does 
not develop over the propagation distance of the experiment. If, on the other hand, 
the instability does develop over such distances, one can say that the soliton is also 
"physically unstable" . 

In order to demonstrate the drift instability associated with violation of the spectral 
condition, and in particular, the importance of the magnitude of Aq , we solve Eq. (J3j) 
with d = 1 and p = 3 for a sinusoidal lattice 

V(Nx) = V cos(2nNx), (54) 
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and also for a Kronig-Penney lattice with the unit cell that consists of a periodic array 
of cells of size 1/N, where for each cell, 




\X\ ^ A A7 



V(Nx) = 7' 1 ^ *«. 1 (55) 



47V < Fl < 2JV 



We excite the instability by shifting the soliton center slightly off the lattice maximum, 
i.e., we use the initial condition A (x) = Uv{x — 8 C ). In Fig. [7J we show the center 
of mass of the solution for N = 0.07, v = 10, Vq = 2.5 and 5 C = 10~ 4 . For these 
parameters, (x(0)) = 5 C and (x(0)) = so that by Eq. f[53l . 



(x,(z)> ~ 5 c cosh(^), = 2\/ ^° \ (56) 

This exponential drift-rate is indeed observed in the simulation for the sinusoidal lattice 
soliton, see Fig. [7J This shows that while the sign of A determines whether the soliton 
is ("mathematically") stable or unstable, the magnitude of | Aq | determines the rate of 
the instability dynamics. 

The drift rate for the KP lattice soliton is several orders of magnitude smaller than 
for the sinusoidal lattice soliton. Intuitively, this is because unlike the sinusoidal lattice, 
the KP lattice affects the soliton profile (and hence the dynamics) only in the soliton 
tail region. As expected, the magnitude of Aq is much larger for the sinusoidal lattice 
soliton {\ { N) = -0.05) than for the KP lattice soliton (A^ ^ -2 ■ 10" 5 ). Moreover, 
the drift rate of the KP lattice soliton is considerably smaller than the one predicted by 
Eq. ( 1561) with A^ = —2 • 10 -5 . This "mismatch" is not surprising, since Eq. ( 1561) is not 
valid for the KP lattice, see also Section 13.11 

At a propagation distance of z — 5, both the sinusoidal and the KP lattice solitons 
hardly shift from their initial location, see Fig. [HJ At a propagation distance of z = 10, 
however, the sinusoidal lattice soliton drifts more than one soliton width whereas the 
Kronig-Penney lattice soliton hardly drifts at all. In that sense, since the propagation 
distance in the simulations corresponds to a distance of 20 diffraction lengths, which 
is longer than most devices in optics, the "mathematically unstable" KP soliton is 
"physically stable". 



6. Discussion and comparison with previous studies 

Most rigorous studies on stability and instability of lattice solitons are based on the 
Grillakis, Shatah and Strauss (GSS) theory [521 [53]. Let Uv > 0, let 



d(u) = H + uV 



|V^)| 2 + (V(Nx lat ) + v) (u^) 2 - H N) Y +1 

' p + 1 



let p(d") = 1 if d" > and p(d") = if d" < 0, and let n_(Z#J) be the number 
of negative eigenvalues of the operator Then, Uv e wz is orbitally stable if 

n_(L^2) = p(d"), and orbitally unstable if n_(L^2) —p(d") is odd [52], [53] . For example, 
stability of lattice solitons was studied in [SH ESI ESI ES] using the GSS theory. In 
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Figure 7. Center of mass of the solution of Eq. ([3]) with d = 1, p = 3 and a sinusoidal 
lattice ([M)) (solid line) and a KP lattice ([55| (dashed line). The lattice parameters 
are N = 0.07 and Vq = 2.5; the initial shift of the soliton center is 6 C = 10 -4 . The 
analytical formula ([5u]) (red dots) is nearly indistinguishable from the numerical result. 



z = z = 5 z = 10 




Figure 8. Beam profiles at several propagation distances for the data of Fig. [7] The 
beam profiles for the sinusoidal lattice (solid line) and the KP lattice (|55[) (dashed 
line) at z — and z — 5 are indistinguishable. 

addition, after this paper was submitted, we found out that the GSS theory was applied 
to narrow lattice solitons in the critical case by Lin and Wei |34j. 

Since d'{v) = J (uv^j dx, the sign of d" is the same as the sign of the power slope. 
Hence, in the GSS theory stability and instability depend on a combination of the slope 
condition and a spectral condition: If both the slope condition and the spectral condition 
are satisfied, the soliton is stable, whereas if either the slope condition is satisfied and 
n_(L^2) is even, or if the slope condition is violated and n_(l/^2) is odd, the soliton is 
unstable. There are two cases not covered by the GSS theory: When the slope condition 
is satisfied and n-(L+l) is odd, and when the slope condition is violated and n-(L+2) 
is even. As Theorem I4.ll shows, in both cases the solitons are unstable. Hence, there is 
a "decoupling" of the slope and spectral conditions, in the sense that both are needed 
for stability, and violation of either of them would lead to instability. 

In [301 EH HO] it was observed numerically that violation of the slope condition 
leads to a width instability, whereas violation of the spectral condition leads to a drift 
instability. Unlike these studies, in this study we prove that violation of the spectral 
condition leads to a drift instability. Moreover, we show that a drift instability occurs 
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in any direction Xj for which the corresponding eigenvalue An „■ is negative, and that the 
drift rate is determined by the magnitude of 1<|. | This further shows that violatiou 
of the spectral condition leads to an instability, regardless of the slope condition and of 
whether n_(L+2) is even or odd. 

In previous studies it was also observed that in the subcritical case, lattice solitons 
centered at a lattice minimum of all widths are stable. In the critical case, it was shown 
that lattice solitons are stable only if they are narrower than a few lattice periods, 
see e.g., [I7J [19]. These results are in agreement with Table Q] in the subcritical and 
critical cases, and imply that our analytical results are valid beyond the regime of 
narrow lattice solitons. In [20, El] it was also shown that in the supercritical case, the 
lattice can stabilize sufficiently wide lattice solitons centered at a lattice minimum but 
cannot stabilize narrow lattice solitons, in agreement with our results. Note, however, 
that unlike most previous works, our results are valid for any dimension d, lattice 
dimension di at and nonlinearity exponent p. 

Another difference from previous studies on linear lattices is that we introduce a 
quantitative approach to the notions of stability and instability. Thus, we show that 
the strength of radial stabilization depends on the magnitude of the slope. Hence, 
in the critical case, the stability of the soliton is "mathematical" but not "physical". 
Similarly, we show that the strength of the transverse instability depends on the value 
of the perturbed zero eigenvalue Aq . Hence, for narrow solitons centered at a lattice 
maximum, the instability is "mathematical" but not necessarily "physical". In such 
cases, the stabilization/destabilization of narrow lattice solitons is highly sensitive to 
the lattice details. This sensitivity becomes smaller as the soliton width increases, and 
is of considerably less importance for 0(1) solitons, which is probably why this feature 
was not observed in previous studies. 
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Appendix A. Proof of Lemma 13.11 

The approach used here is similar to [301 EH] • Substituting the expansion ffl9]) in Eq. ffl5l) 
gives 

V%v + 4 - (l + A^x^))^ + 0(A 4 ) = 0. (A.l) 

jt A generalization of these results to non-narrow beams can be found in |57j . 
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Let Ujy-(x) be given by Eq. (|T7|) . Then, the equation for g is 
V 2 ^(x) + P U p - 1 g -vg = V 2 (± lat )U . 

Therefore, 

^(x) = -L; 1 [V r 2 (x^)W(|x|)]. (A.2) 
Appendix B. Proof of Lemma 14.11 

By Eq. fl2M . the power of the rescaled lattice soliton = J (uff(x)) 2 dx is given by 
Vf, = V v=x - 2N 2 J U{f)L^[V 2 {x lat )U]d± + 0{N A ) 

= V„=i - 2N 2 J VzfaMr)^ 1 [U] rfx + 0(A> 4 ), (B.l) 

where V v= \ = j U 2 (r)d5t and r = |x|. In order to proceed, we prove the following 
Lemma: 



Lemma Appendix B.l Let U v be the solution of Eq. (27\ ) and let L +tV be given by 
Eq. g§. Then, L^Ur, = -d v U v . 

Proof: Differentiating Eq. ff27j) with respect to i] gives 

d v (V 2 U V ) + d v U* - d v ( V U V ) = V 2 {d n U v ) + pUP- 1 (d v U v ) -U v - V d v U v ■. 



-L + d v U n -U v = 0. □. 



Since U v (f) = rjp- 1 U(- s /fjr), then 

1 _l x _j_ (\ _i \ 

dM„ = tV p1 U + V' 1 -kV 2r ) Uf- 

p-1 \2 J 

Therefore, L+ l U = L^U r]=x = -{d v U v ) v= i = — - \rU f . Substituting in Eq. (IPO) 



gives 



V$ = V v=1 + 2N 2 j V 2 {± lat )U {-^ + f^j dx + 0(N A ). (B.2) 
Since V 2 is given by Eq. ([2Tj) . Eq. (1B.2|) can be written as 

d-lat 

= V u=1 - C V N 2 % + 0(N A ), 



3=1 

where Cy is given by 



C v = - I x 2 ( ^— + fUUf ) d± (B.3) 



p — 1 

To bring Cy to the form ( 1361) . we note that 

V ■ (6(f)x) = -^~i^z (r d b(r)) = ^ {df d ~ l b + f d b') = db + fb'. 
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Substituting 6(f) = f p 31 " W 2 (f) shows that 



V ■ (fi^- d W 2 (r)xJ = dr^ d U 2 + ~ d ) >~~"W - 2r~~" + W 

Thus, we can rewrite Eq. (IB. 3D as 

°y= - \l ri v • f^" d w 2 x) dx (B.4) 
= I/ r A-V X .v(^) rfS 

2 J \fT->- d VP-1 Jf^- d+1 J 

Finally, by the dilation transformation (1141) . 

Vi N) = J (4 iV )(x)) 2 rfx = ^ y (u f[ (±)) 2 d^ = V 4 -W 1 V^. 

Appendix C. Proof of Lemma 14.21 

Consider the eigenvalue problem 

Multiplying Eq. flC.lj) by fjff and integrating gives 

We recall that in the absence of a lattice, the operator L+l reduces to L +jV , see 
Eq. (|42l) . which has d zero eigenvalues Aoj = 0, with the corresponding eigenfunctions 
for j = l,...,cf, see Eq. ( 1431) . By Eq. (|26l) . in the presence of the lattice, 



4 j\ ^^-^7/2 , r,~-^ T -d+l 7n ,, 



u { v N) =U V + r]^0(N 2 ). Similarly, by Eq. (jlBjl. Eq. (fT9l) and Eq. we can expand 
the potential as 

7(iVx /at ) = 7(A>x /at ) = 7(0) + vV{N5c lat ) = 7(0) + V (N 2 V 2 (± lat ) + 0(iV 4 )) (C.3) 
= 7(0) + AT 2 Vjj&j + *7 • C(^ 4 ) = ^(0) + r/iV 2 Vjjtf + V ■ 0(N A ). 

3=1 j=l 

Consequently, the operator L+2 can be expanded as 

L { + N l = - V 2 - p {u^y- 1 + v + V(Nx lat ) (C.4) 
= - V 2 - p (U n + r]^0(N 2 )Y 1 + v + 7(0) + 0{N 2 ) 
= - V 2 - pW^- 1 + 77 + 0(N 2 ) = L +iV + 0(N 2 ). 
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Therefore, we expand 



C(x) = (1 + 0(iV 2 )) , A$ = 5,iV 2 + 0(iV 4 ). (C.5) 

F (JV) 



By Eqs. (I26p and (10.5p . we can also rewrite the eigenfunction y v - as 

/ ^ ,(x) = ^-( 1 + 0(iv ' 2) '- (C ' 6) 

We now use the approximations (1C.5|) and (1C.6|) in order to evaluate the terms in 
Eq (|C3j» . By Eq. (1C3D . the right-hand-side of Eq. flU^j) is equal to 

Ag } / (/r } ) 2 rfx= (iV 2 5, + 0(N*)) (J (^) 2 ^ + 0(iV 2 

/ (|^x+o(in (c.7) 

By Eq. ( 1C.6I) the left-hand-side of Eq. ( 1C2[) . approximation ( 1C6f) is equal to 

/ f%>L™f™dx = J ^-L^^-dx + 0(iV 4 ), (C.8) 



r4\ 



where the error term is 0{N ) due to the properties of the Rayleigh quotient, see 
e.g., [SB]. 

The integral term on the right-hand-side of Eq. (1C.8|) is equal to 

/¥ L "C dx ^/K K, ) 2 ^(iVx,Mx. (C) 

Indeed, differentiating Eq. (TT3T) with respect to iCj gives 

L W^ = _/ ^(iVx to ) \ w (ai0) 

Multiplying Eq. fIC.lOl) by ^-ui N \ integrating over x and integrating by parts gives 
Eq. flCT9l) . Using Eq. fl03l . the right-hand-side of Eq. fl09l) is given by 

\ J H N) ) 2 ^2 V ( N ^t)d^ = V N 2 v n J Ufa+0(N 4 ). (C.ll) 

Comparing the approximation flC.7j) for the left-hand-side of Eq. (1C.2|) with the 
approximation ( IC.llj) for the right-hand-side of Eq. ( 1C.2|) shows that 



s ' I = J (C12) 

5 i = Wjj /ni \ 2 = dv 3j jj^ . (C. 13) 



Hence, 



r f 9U A' 

J \dx 3 ) 



Similar results were obtained in [31] for a soliton centered at a general non-degenerate 
critical point of the lattice (i.e., without assuming that the critical point is symmetric 
with respect to x[°]). 
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By the Pohozaev identities for Eq. (fl8l) (see [59J, pp. 76) 
Therefore, we get that 

Sj = Svjj. 



fu 2 



p(2-ri)+2+rf 
d(p-l) 



_ S 



— d' 



(C.14) 



Appendix D. Computing small eigenvalues of a very large matrix 

When d > 2, the discretized operator is represented by an extremely large matrix. 
Hence, straightforward application of standard numerical routines (such as Matlab's 
eig/eigs) usually either fails to give accurate results or does not converge. 

In order to overcome this numerical problem, we used a more efficient and robust 
numerical method based on the Arnoldi algorithm (performed by ARPACK [60] . which 
is available in Matlab through the function eigs). Essentially, we compute the largest- 
magnitude eigenvalues of the inverse matrix A^ 1 which correspond to the smallest 
eigenvalues of the matrix A. 

We compute the LU factorization of A with complete pivoting. Then, we shift the 
values on the main diagonal of U by a small value in order to avoid numerical errors that 
might result from singularity of the matrix during the computation of A~ l . Then, in 
order to avoid working with the explicit from of the inverse matrix A^ 1 which is dense, 
we compute A -1 implicitly through the subfunction LUPinv and apply it to the function 
eigs. This way, we exploit the sparsity of the LU factorized matrices U and L. The 
function eigs then computes the desired number of eigenvalues of largest magnitude. 

The following code was given to us by Prof. S. Toledo: 

function [V,d] = ev_calculation(A,ev_number,eps) 
[m n] = size(A); normA = norm(A,l); 
[L.U.P.Q] = lu(A.l.O); 
for j=l:n 

if (abs(U(j , j)) < eps*normA) 

U( j , j ) = eps*normA; 
end 
end 

h = @LUPinv; 
opts.issym = true; 
opts.isreal = true; 
opts.tol = eps; 

[V,D] = eigs(h,n,ev_number, 'LM ; ,opts) ; 



function Y = LUPinv (X) 
Yl = P*X; 
Y2 = L \ Yl; 
Y3 = U \ Y2; 
Y = Q*Y3; 
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end 

end 



Appendix E. Proof of Lemma 14.31 

Multiplying Eq. (jHJ) by A* and subtracting the conjugate equation gives 

^-\A\ 2 = tA*V 2 A + c.c, (E.l) 
dz 

where c.c. stands for complex conjugate. Multiplying by x and integrating over x gives 
^ J x|A| 2 = J ixA*V 2 A + c.c. = -% J VA(dA* + x ■ VA*) + c.c. 

= 2dlm J A*VA. (E.2) 
Differentiating Eq. ( IE. 21) yields 

J x|A| 2 = 2d Im j (A* Z VA + A*VA Z ) 

= 2dlm J (A*^ A - A Z VA*) = Ad Im J A* Z VA 

= - Ad Re J (V 2 A* + lA^A* - V{Nx)A*) VA. 

The first two terms vanish since they are complete derivatives. Therefore, 
d 2 



dz 



J x|A| 2 = 4dRe J V(Nx)A*VA 

= 2d J V(Nx lat )V\A\ 2 = -2d J \A\ 2 VV(Nx lat ). (E.3) 



Finally, by Eq. (J22]) 

dz 2 
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